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Abstract. In semiclassical theories for chaotic systems such as Gutzwiller's periodic 
orbit theory the energy eigenvalues and resonances are obtained as poles of a non- 
convergent series g(w) = A n exp(is„w). We present a general method for the 
analytic continuation of such a non-convergent series by harmonic inversion of the 
"time" signal, which is the Fourier transform of g(w). We demonstrate the general 
applicability and accuracy of the method on two different systems with completely 
different properties: the Riemann zeta function and the three disk scattering system. 
The Riemann zeta function serves as a mathematical model for a bound system. We 
demonstrate that the method of harmonic inversion by filter-diagonalization yields 
several thousand zeros of the zeta function to about 12 digit precision as eigenvalues of 
small matrices. However, the method is not restricted to bound and ergodic systems, 
and does not require the knowledge of the mean staircase function, i.e., the Weyl 
term in dynamical systems, which is a prerequisite in many semiclassical quantization 
conditions. It can therefore be applied to open systems as well. We demonstrate this 
on the three disk scattering system, as a physical example. The general applicability 
of the method is emphasized by the fact that one does not have to resort a symbolic 
dynamics, which is, in turn, the basic requirement for the application of cycle expansion 
techniques. 
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1. Introduction 

Since the development of periodic orbit theory by Gutzwiller |], @ it has become a 
fundamental question as to how individual semiclassical eigenenergies and resonances 
can be obtained from periodic orbit quantization for classically chaotic systems. A major 
problem is the exponential proliferation of the number of periodic orbits with increasing 
period, resulting in a divergence of Gutzwiller's trace formula at real energies and below 
the real axis, where the poles of the Green's function are located. The periodic orbit 
sum is a Dirichlet series 

g(w) = J2 A nj SnW , (1) 

n 

where the parameters A n and s n are the amplitudes and periods (actions) of the periodic 
orbit contributions. In most applications Eq. [I] is absolutely convergent only in the region 
Im w > c > with c the entropy barrier of the system, while the poles of g(w), i.e., the 
bound states and resonances, are located on and below the real axis, Imw < 0. Thus, to 
extract individual eigenstates, the semiclassical trace formula (|lj) has to be analytically 
continued to the region of the quantum poles. 

Up to now no general procedure is known for the analytic continuation of a non- 
convergent Dirichlet series of the type of Eq. [I]. All existing techniques are restricted to 
special situations. For bound and ergodic systems the semiclassical eigenenergies can 
be extracted with the help of a functional equation and the mean staircase function 
(Weyl term), resulting in a Riemann-Siegel look-alike formula 0, ^, |^, Efl. Alternative 
semiclassical quantization conditions based on a semiclassical representation of the 
spectral staircase 0, [5] and derived from a quantum version of a classical Poincare 
map H are also restricted to bound and ergodic systems. 

For systems with a symbolic dynamics the periodic orbit sum (P can be reformulated 
as an infinite Euler product, which can be expanded in terms of the cycle length of the 
symbolic code. If the contributions of longer orbits are shadowed by the contributions 
of short orbits the cycle expansion technique can remarkably improve the convergence 
properties of the series and allows to extract the bound states and resonances of bound 



and open systems, respectively [10, 11, 12, 13 1 . A combination of the cycle expansion 



technique with a functional equation for bound systems has been studied by Tanner 



et al. jlj]. However, the existence of a simple symbolic code is restricted to very few 
systems, and cycle expansion techniques cannot be applied, e.g., to the general class of 
systems with mixed regular-chaotic classical dynamics. 

In this paper we present a general technique for the analytic continuation and the 
extraction of poles of a non- convergent series of the type of Eq. [I]. The method is 
based on harmonic inversion by filter-diagonalization. The advantage of the method 
is that it does not depend on special properties of the system such as ergodicity or 
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the existence of a symbolic dynamics for periodic orbits. It does not even require the 
knowledge of the mean staircase function, i.e., the Weyl term in dynamical systems. 
The only assumption we have to make is that the analytic continuation of the Dirichlet 
series g(w) (Eq. [I]) is a linear combination of poles (w — u^) -1 , which is exactly the 
functional form of, e.g., a quantum mechanical response function with real and complex 
parameters Wk representing the bound states and resonances of the system, respectively. 
To demonstrate the general applicability and accuracy of our method we will apply it 
to two systems with completely different properties, first the zeros of the Riemann zeta 



function |y| [16[], as a mathematical model for a bound system, and second the three 
disk scattering system as a physical example. 

As pointed out by Berry || the density of zeros of Riemann's zeta function can 
be written, in formal analogy with Gutzwiller's semiclassical trace formula, as a non- 
convergent series, where the "periodic orbits" are the prime numbers. A special property 
of this system is the existence of a functional equation which allows the calculation of 
Riemann zeros via the Riemann-Siegel formula [R| [L7|]. An analogous functional 
equation for quantum systems with an underlying chaotic (ergodic) classical dynamics 
has served as the basis for the development of a semiclassical quantization rule for bound 
ergodic systems H f|, |^, |6| . The Riemann zeta function has also served as a mathematical 
model to study the statistical properties of level distributions jI7|, [L8|, [HJ. We will 
demonstrate that harmonic inversion can reveal the Riemann zeros with extremely 
high accuracy and with just prime numbers as input data. The most important 
advantage of our method is, however, its wide applicability, i.e., it can be generalized in 
a straightforward way to non-ergodic bound or open systems. 

Our second example, the three disk scattering problem, is an open and non-ergodic 
system. Its classical dynamics is purely hyperbolic, and the periodic orbits can be 
classified by a complete binary symbolic code. This system has served as a model 
for the development of cycle expansion techniques JlO], 0, |13|| . When applying the 



harmonic inversion technique to the three disk scattering system we will highlight the 
general applicability of our method by not having to make use of its symbolic dynamics 
in any way. 

It is evident that methods invoking special properties of a given system may be more 
efficient regarding, e.g., the number of periodic orbits required for the calculation of a 
certain number of poles of the response function g(w) in that particular case. It is not our 
purpose to compete with the efficiency of such methods. Rather, the advantage of the 
harmonic inversion technique lies in its wide applicability, which allows the investigation 
also of systems not possessing special properties. This is demonstrated in this paper by 
solving two completely different problems, viz. the zeros of the Riemann zeta function 
and the three disk scattering system, with one and the same method. 
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The paper is organized as follows. In Section 2 we explain the general idea of the 
method by way of example of the Riemann zeros. This is followed by the derivation of 
the harmonic inversion method in Section 3 and the presentation of numerical results for 
the Riemann zeros in Section 4. The method is extended to the general case of periodic 
orbit quantization in Section 5, and its usefulness and wide applicability is demonstrated 
for the three disk scattering system, as a physical example, in Section 6. 

2. The Riemann zeta function 

Our goal is to introduce our method for periodic orbit quantization by harmonic 
inversion using, as an example, the well defined problem of calculating zeros of the 
Riemann zeta function. There are essentially two advantages of studying the zeta 
function instead of a "real" physical bound system. First, the Riemann analogue of 
Gutzwiller's trace formula is exact, as is the case for systems with constant negative 
curvature || [8] , whereas the semiclassical trace formula for systems with plane geometry 
is correct only to first order in h. This allows a direct check on the precision of the 
method. Second, no extensive periodic orbit search is necessary for the calculation of 
Riemann zeros, as the only input data are just prime numbers. It is not our intention 
to introduce yet another method for computing Riemann zeros, which, as an objective 
in its own right, can be accomplished more efficiently by specific procedures. Rather, 
in our context the Riemann zeta function serves primarily as a mathematical model to 
illustrate the power of our technique when applied to bound systems. 



2.1. General remarks 

Before discussing the harmonic inversion method we start with recapitulating a few 
brief remarks on Riemann's zeta function necessary for our purposes. The hypothesis 
of Riemann is that all the non-trivial zeros of the analytic continuation of the function 

((z) = n ~ Z = II ( 1 ~ P Z ) ' ( Re z > P : P rimes ) ( 2 ) 

n=l p 

have real part |, so that the values w = Wk, defined by 

C Q - iw^ = 0, (3) 

are all real or purely imaginary Jl5, 16| . The Riemann staircase function for the zeros 
along the line z = ~ — iw, defined as 

oo 

N{w) = y £ t Q{w-w k ), (4) 

k=l 



5 



i.e. the number N(w) of zeros with < w, can be split [||, [1^, [IIJ into a smooth part, 
N(w) = — arg T ( \ + \%w \ - ^- In 7r + 1 

7T \4 2 / 27T 

= ™ flnj^l - lU - + ^ — % + 0( W S) , (5) 

and a fluctuating part, 

N osc (w) = -- lim Im ln( (- - i{w + 17])] . (6) 

Substituting the product formula (||) (assuming that it can be used when Re z = |) 
into (Q) and expanding the logarithms yields 

1 OO 1 

N osc (w) = -- Im E E — ^ ^ mln(p) • (7) 

" p m=l l,L t ) 

Therefore the density of zeros along the line z = | — iw can formally be written as 

£oscM = = —- lmg(w) (8) 

with the response function g(w) given by the series 

p m=l V 

which converges only for Im w > | . 

Obviously Eq. § is of the same type as the response function ([[]), with the entropy 
barrier c = ~, i.e., Eq. |] does not converge on the real axis, where the Riemann zeros 
are located. The mathematical analogy between the above equation and Gutzwiller's 
periodic orbit sum 

g OBC (E) « Im £ „4 po e iS »°, (10) 

po 

with Ap the amplitudes and S^o the classical actions (including phase information) 
of the periodic orbit contributions, was already pointed out by Berry j3], |J. For the 
Riemann zeta function the primitive periodic orbits have to be identified with the primes 
p, and the integer m formally counts the "repetitions" of orbits. The "amplitudes" and 
"actions" are then given by 

A =^ (11) 

•^P m pm/2 ' ^ > 

S pm = mw\n(p) . (12) 

Both equation (|8]) for the Riemann zeros and - for most classically chaotic physical 
systems - the periodic orbit sum ( PTJ|) do not converge. In particular, zeros of the 
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zeta function, or semiclassical eigenstates, cannot be obtained directly using these 
expressions. The problem is to find the analytic continuation of these equations to 
the region where the Riemann zeros or, for physical systems, the eigenenergies and 
resonances, are located. Eq. |9] is the starting point for our introduction and discussion 
of the harmonic inversion technique for the example of the Riemann zeta function. The 
generalization of the method to periodic orbit quantization (Eq. |1(|) in Section 5 will be 
straightforward. 

Although Eq. ^| is the starting point for the harmonic inversion method, for 
completeness we quote the Riemann-Siegel formula, which is the most efficient approach 
to computing Riemann zeros. For the Riemann zeta function it follows from a functional 
equation [15| that the function 



Z(w) = exp < —i 



,1 1. \ 1 ' 

argr — I — iw wmn 

1 4 2 J 2 



}c(i-*.) (13) 



is real, and even for real w. The asymptotic representation of Z(w) for large w, 



, . _ ^ cos{tcN(w) - wlnn} 
Z(w)--2 ^ — - 2 



n=l " 
1 



_ r^tlV^) ( 2 _1) i cos( 27 r(t 2 -t-l/16)) + _ _ _ 
V w J cos (2irt) 

with t = \Jw /2ix — Int[y^y/2~7r] is known as the Riemann-Siegel formula and has been 
employed (with several more correction terms) in effective methods for computing 
Riemann zero s ||17||. Note that the principal sum in (14) has discontinuities at integer 



positions of ^w/2ir, and therefore the Riemann zeros obtained from the principal sum 
are correct only to about 1 to 15 percent of the mean spacing between the zeros. The 
higher order corrections to the principal Riemann-Siegel sum remove, one by one, the 
discontinuities in successive derivatives of Z{w) at the truncation points and are thus 
essential to obtaining accurate numerical results. An alternative method for improving 
the asymptotic representation of Z(w) by smoothing the cut-offs with an error function 
and adding higher order correction terms is presented in 0. 

An analogue of the functional equation for bound and ergodic dynamical systems 
has been used as the starting point to develop a "rule for quantizing chaos" via a 
"Riemann-Siegel look-alike formula" 0, |5], H|. This method is very efficient as it requires 
the least number of periodic orbits, but unfortunately it is restricted to ergodic systems 
on principle reasons, and cannot be generalized either to systems with regular or mixed 
classical dynamics or to open systems. 

By contrast, the method of harmonic inversion does not have these restrictions. We 
will demonstrate that Riemann zeros can be obtained directly from the "ingredients" 
of the non- convergent response function (^), i.e., the set of values A pm and S pm , thus 
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avoiding the use of the functional equation, the Riemann-Siegel formula, the mean 
staircase function (5), or any other special property of the zeta function. The comparison 
of results in Section 4 will show that the accuracy of our method goes far beyond the 
Riemann-Siegel formula (14) without higher order correction terms. The main goal of 
this paper is to demonstrate that because of the formal equivalence between Eqs. (^) 
and ( |T0D our method can then be applied to periodic orbit quantization of dynamical 
systems [JZIJ without any modification. 



2.2. The ansatz for the Riemann zeros 

To find the analytic continuation of Eq. (|9|) in the region Im w < \ we essentially wish 
to fit g(w) to its exact functional form, 

aUw) = — dk , ■ > (15) 

arising from the definition of the Riemann staircase (01). The "multiplicities" dk in Eq. 



C(s) = — g(w)e-^dw = ij:j:-^8(s-mln(p)), (16) 

Z7T J —oo „ ™_i V 



15 are formally fitting parameters, which here all should be equal to 1. 

It is hard to directly adjust the non- convergent (on the real axis) series g(w) to the 
form of Qex(w). The first step towards the solution of the problem is to carry out the 
adjustment for the Fourier components of the response function, 

p m=l f 

which after certain regularizations (see below) is a well-behaved function of s. Due to 
the formal analogy with the results of periodic orbit theory (see Eqs. 11 and 12), C(s) 
can be interpreted as the recurrence function for the Riemann zeta function, with the 
recurrence positions S pm = m ln(p) and recurrence strengths of periodic orbit returns 
Apm = iln(p)p~ m / 2 . The exact functional form which now should be used to adjust 
C(s) is given by 



1 r+OO 00 

C cx (s) = — / g ex (w)e- isw dw = -i £ d k e~ iw * s . (17) 
2tt J-co k=1 



oo 



C cx (s) is a superposition of sinusoidal functions with frequencies f Wk given by the 
Riemann zeros and amplitudes dk — 1. 

Fitting a signal C (s) to the functional form of Eq. ([H]) with, in general, both complex 
frequencies Wk and amplitudes dk is known as harmonic inversion, with a large variety 



of applications in various fields [^T|, ^3], 24, |25|]. The harmonic inversion analysis is 



especially non-trivial if the number of frequencies in the signal C(s) is large, e.g., more 

f It is convenient to use the word "frequencies" for Wk referring to the sinusoidal form of C(s). We 
will also use the word "poles" in the context of the response function g(w). 
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than a thousand. It is additionally complicated by the fact that the conventional way 
to perform the spectral analysis by studying the Fourier spectrum of C(s) will bring us 
back to analyzing the non- convergent response function g(w) defined in Eq. [| Until 
recently the known techniques of spectral analysis |^TJ would not be applicable in the 
present case. It is the filter-diagonalization method [22, |2J| [2J1] which have turned the 
harmonic inversion concept into a general and powerful computational tool. 

The signal C(s) as defined by Eq. [IB] is not yet suitable for the spectral analysis. 
The next step is to regularize C(s) by convoluting it with a Gaussian function to obtain 
the smoothed signal, 

C a (s) = -^L- / +0 ° C(s')e~^ 2 ^ 2 ds' 
v27rcr J-c 



oo 



EE 



Xl \n(p) 



27R7 V^l^ 2 



e 



-(s-mln(p)) 2 /2a 2 



that has to be adjusted to the functional form of the corresponding convolution of C ex (s). 
The latter is readily obtained by substituting d k in Eq. [17| by the damped amplitudes, 

d k ^d[ a) = d k e- w ' u2/2 . (19) 

The regularization (18) can also be interpreted as a cut of an infinite number of high 
frequencies in the signal which is of fundamental importance for numerically stable 
harmonic inversion. Note that the convolution with the Gaussian function is no 
approximation, and the obtained frequencies w k and amplitudes d k corrected by Eq. 



19| are still exact, i.e., do not depend on a. The convolution is therefore not related to 
the Gaussian smoothing devised for Riemann zeros in [2f| and for quantum mechanics 
in [27j, which provides low resolution spectra only. 

Before proceeding further we note that even though the derivation of Eq. 18 assumed 
that the zeros w k are on the real axis, the analytic properties of C a (s) imply that its 
representation by Eq. 18 includes not only the non-trivial real zeros, but also all the 
trivial ones, w k = —i(2k + |), k — 1,2, . . ., which are purely imaginary. The general 
harmonic inversion procedure described below does not require the frequencies to be 
real. Both the real and imaginary zeros w k will be obtained as the eigenvalues of a 
non-Hermitian generalized eigenvalue problem. 



3. Filter-diagonalization method for harmonic inversion 

The harmonic inversion problem can be formulated as a non-linear fit (see, e.g., Ref. 
nj) of the signal C(s) defined on an equidistant grid, 

c n = C(nr)=J2d k e- inTWk , n = 0, 1, 2, . . . N, (20) 
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with the set of generally complex variational parameters {wk, (In this context the 
Discrete Fourier Transform scheme would correspond to a linear fit with N amplitudes 
dk and fixed real frequencies Wk = 2nk/Nr, k = 1,2, . . .N. The latter implies the so 
called "uncertainty principle", i.e., the resolution, defined by the Fourier grid spacing, 
Aw, is inversely proportional to the length, s max = Nt, of the signal C(s).) The "high 
resolution" property associated with Eq. |20| is due to the fact that there is no restriction 
for the closeness of the frequencies W}. as they are variational parameters. In Ref. 



it was shown how this non-linear fitting problem can be recast as a linear algebraic one 
using the filter-diagonalization procedure. The essential idea is to associate the signal 
c n with an autocorrelation function of a suitable dynamical system, 



(21) 



where ( • , • ) defines a complex symmetric inner product (i.e., no complex conjugation). 
The evolution operator can be defined implicitly by 



K 



U = e~ iTQ = J2e~ iT " k \Tk)(Tk\ , (22) 
k=i 

where the set of eigenvectors {T^} is associated with an arbitrary orthonormal basis set 
and the eigenvalues of U are = e~ %TUJk (or equivalently the eigenvalues of Cl are u^). 
Inserting Eq. |22| into Eq. |21] we obtain Eq. with 

4 = (T fc ,$ ) 2 , (23) 

which also implicitly defines the "initial state" $o- 

This construction establishes an equivalence between the problem of extracting 
spectral information from the signal with the one of diagonalizing the evolution operator 
U = e~ lrn (or the Hamiltonian Cl) of the fictitious underlying dynamical system. 
The filter-diagonalization method is then used for extracting the eigenvalues of the 
Hamiltonian Cl in any chosen small energy window. Operationally this is done by 
solving a small generalized eigenvalue problem whose eigenvalues yield the frequencies 
in a chosen window. The knowledge of the operator Cl itself is not required, as for a 
properly chosen basis the matrix elements of Cl can be expressed only in terms of c n . The 
advantage of the filter-diagonalization procedure is its numerical stability with respect to 
both the length and complexity (the number and density of the contributing frequencies) 
of the signal. Here we apply the method of Ref. |23[] which is an improvement of the 
filter-diagonalization method of Ref. |22| in that it allows to significantly reduce the 



required length of the signal by implementing a different Fourier-type basis with an 
efficient rectangular filter. Such a basis is defined by choosing a small set of values 
ifj in the frequency interval of interest, Tu} min < <pj < ru max , j = 1,2,..., J, and the 
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maximum order, M, of the Krylov vectors, $ n = e m ™ 2 <|> , used in the Fourier series, 

M M 

m j = m( Vj ) = Y J e intp ^ n = e m{ ^- TQ) %. (24) 

n=0 n=0 

It is convenient to introduce the notations, 

U$ = U^(^, Vr ) = (*(<^), e"^%(^)) , (25) 

for the matrix elements of the operator e~ iprn , and XJ^ P \ for the corresponding small 
J x J complex symmetric matrix. As such denotes the matrix representation of the 
operator U itself and U^ -*, the overlap matrix with elements (^(tpj), ty((pj>)), which is 
required as the vectors ^f((pj) are not generally orthonormal. Now using these definitions 
we can set up a generalized eigenvalue problem, 

U (p) B fc = e^ p ™ fc U (0) B fc , (26) 

for the eigenvalues e~ tpTUJk of the operator e~ ipT . The column vectors with elements 
Bjk define the eigenvectors in terms of the basis functions tyj as 

j 

T k = J2 B ^j, (27) 

3=1 

assuming that the ^/s form a locally complete basis. 

The matrix elements ( p5[) can be expressed in terms of the signal c n , the explicit 
knowledge of the auxiliary objects Cl, or $ is n °t needed. Indeed, insertion of Eq. 
24| into Eq. use of the symmetry property, (^/, U&j — (ll^, and the definition 
of c n , Eq. ETL gives after some arithmetics 

M 

UM(<p, cp') = {e-* - e-^'y 1 [e-* ]T e m *' c n+p (28) 

ra=0 

M 2M 

_ e -V e irvp c n+p - e iMip e i( - n ~ M ~^ v ' 

n=0 n=M+l 
2M 

n=M+l 

2M 

' ft it I 7i /r „ I i 1 \ „inifi„ 

n+p- 



U®((p,(p) = £(M- \M-n\ + l)e^c 



n=0 

(Note that the evaluation of requires knowledge of c n for n — p,p + 1, . . . , N — 
2M + p.) 

The solution of the generalized eigenvalue problem fl26|) is usually done by a singular 
value decomposition of the matrix U^ -*. Each value of p yields a set of frequencies Wk 
and, due to Eqs. [23|, |24| and ^7|, amplitudes, 

4= (j2 B jkflc n e in A . (29) 

Vi=l n=0 J 



11 



Note that Eq. |29] is a functional of the half signal c n , n = 1,2,..., M. Even though 
in all our applications Eq. [2^ yields very good results, here we present an even better 
expression for the coefficients dk (see Ref. p4|), 

2 



di 



M + 
1 

W- 



1 J 



j 



-Y.B jk U^\^,Wk) 

1 7=1 



(30) 



with U^((p j} w k ) defined by Eq. 28. Eq. 30 is a functional of the whole available signal 



c n , n = 0,l,...,2M. 



The converged Wk and dk should not depend on p. This condition allows us to 
identify spurious or non-converged frequencies by comparing the results with different 
values of p (e.g., with p = 1 and p = 2). We can define the simplest error estimate e 
as the difference between the frequencies Wk obtained from diagonalizations with p = 1 
and p = 2, i.e. 



\w 



(p=i) 



(P=2), 



(31) 



4. Riemann zeros by harmonic inversion 

We have introduced in Section 2 the method for periodic orbit quantization by harmonic 
inversion for the example of the Riemann zeta function because this model allows a 
direct check of the precision of our method. In this Section we present and discuss 
the numerical results obtained for the Riemann zeros. We also discuss the amount of 
periodic orbit input data, viz. here the prime numbers, required to obtain converged 
semiclassical eigenenergies and Riemann zeros, respectively. 



4-1. Numerical results 

For a numerical demonstration we construct the signal C a (s) using Eqs. [16| and 18 in 
the region s < ln(10 6 ) = 13.82 from the first 78498 prime numbers and with a Gaussian 
smoothing width o = 0.0003. Parts of the signal are presented in Fig. 1. Up to s ~ 8 
the Gaussian approximations to the ^-functions do essentially not overlap (see Fig. 
la) whereas for s ^> 8 the mean spacing As between successive 5-functions becomes 
much less than the Gaussian width a = 0.0003 and the signal fluctuates around the 
mean C(s) = ie s ^ 2 (see Fig. lb). From this signal we were able to calculate about 
2600 Riemann zeros to at least 12 digit precision. For the small generalized eigenvalue 
problem ( |26| ) we used matrices with dimension J < 100. Some Riemann zeros Wk, the 
corresponding amplitudes dk, and the estimated errors e (see Eq. Q3"T|)) are given in 
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Table 1. Within the numerical error the Riemann zeros are real and the amplitudes are 
consistent with d k = 1 for non-degenerate zeros. To fully appreciate the accuracy of 
our harmonic inversion technique we note that zeros obtained from the principal sum of 
the Riemann- Siegel formula (14) deviate by about 1 to 15 percent of the mean spacing 
from the exact zeros. Including the first correction term in (14) the approximations 
to the first five zeros read w 1 = 14.137, w 2 = 21.024, w 3 = 25.018, w A = 30.428, 
and u> 5 = 32.933, which still significantly deviates from the exact values (see Table 1). 
Considering even higher order correction terms the results will certainly converge to 
the exact zeros. However, the generalization of such higher order corrections to ergodic 
dynamical systems is a nontrivial task and requires, e.g., the knowledge of the terms in 
the Weyl series, i.e., the mean staircase function after the constant || 29]. The perfect 
agreement of our results for the w k with the exact Riemann zeros to full numerical 
precision is remarkable and clearly demonstrates that harmonic inversion by filter- 
diagonalization is a very powerful and accurate technique for the analytic continuation 
and the extraction of poles of a non-convergent series such as Eq. [1]. 

A few Wk have been obtained (see Table 2) which are definitely not located on the 
real axis. Except for the first at w = i/2 they can be identified with the trivial real zeros 
of the zeta function at z — —2n\ n — 1, 2, ... In contrast to the nontrivial zeros with real 
Wk, the numerical accuracy for the trivial zeros decreases rapidly with increasing n. The 
trivial zeros w n = —i(2n + |) are the analogue of resonances in open physical systems 
with widths increasing with n. The fact that the trivial Riemann zeros are obtained 
emphasizes the general applicability of our method and demonstrates that periodic orbit 
quantization by harmonic inversion can be applied not only to closed but to open systems 
as well. The decrease of the numerical accuracy for very broad resonances is a natural 



numerical consequence of the harmonic inversion procedure |23|, |24 



The value w = i/2 in Table 2 is special because in this case the amplitude is negative, 



i.e., dk = —1. Writing the zeta function in the form |T6 



C{\ -iw) = C ]J(w - w k ) d *A(w, w k ) (32) 
1 k 

where C is a constant and A a regularizing function which ensures convergence of the 
product, integer values d k are the multiplicities of zeros. Therefore it is reasonable to 
relate negative integer values with the multiplicities of poles. In fact, ((z) has a simple 
pole at z = \ — iw = 1 consistent with w = i/2 in Table 2. 



4-2. Required signal length 

We have calculated Riemann zeros by harmonic inversion of the signal C a (s) (Eq. 18) 
which uses prime numbers as input. The question arises what are the requirements on 
the signal C a (s), in particular what is the required signal length. In other words, how 
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many Riemann zeros (or semiclassical eigenenergies) can be converged for a given set of 
prime numbers (or periodic orbits, respectively). The answer can be directly obtained 
from the requirements on the harmonic inversion technique. In general, the required 
signal length s max for harmonic inversion is related to the average density of frequencies 
g(w) by @ 

■Vax ~ 4:1Tq(w) . (33) 



From Eq. 33 the required number of primes (or periodic orbits) can be directly estimated 
as {# primes p | hip < s max } or periodic orbits | s po < s max }. For the special 
example of the Riemann zeta function the required number of primes to have a given 
number of Riemann zeros converged can be estimated analytically. With the average 
density of Riemann zeros derived from (5), 

g( w ) = — = — In — (34) 



dw 2n \2n / 

we obtain 

2 

.271/ \ Z7T, 

The number of primes with p < p max can be estimated from the prime number theorem 



ln(p m a X ) = 2 In ( ^- J => p max = (^-) ■ (35) 



, v Pmax (W /27T 2 

7T Pmax ~ T~ ( T = . , , v • 36 

ln(p m ax) 2 ln(w/27r) 

On the other hand the number of Riemann zeros as a function of w is given by Eq. (5). 
The estimated number of Riemann zeros which can be obtained by harmonic inversion 
from a given set of primes is presented in Fig. 2. For example, about 80 zeros (w < 200) 
can be extracted from the short signal C a (s) with s max = ln(1000) = 6.91 (168 prime 
numbers) in agreement with the estimates given above. Obviously, in the special case 
of the Riemann zeta function the efficiency of our method cannot compete with that 
of the Riemann-Siegel formula method (14) where the number of terms is given by 
^max = Int [^Jw/2tt] and, e.g., 5 terms in Eq. 14 would be sufficient to calculate good 
approximations to the Riemann zeros in the region w < 200. Our primary intention 
is to introduce harmonic inversion by way of example of the zeros of the Riemann 
zeta function as a general tool for periodic orbit quantization, and not to use it as 
an alternative method for solving the problem of finding most efficiently zeros of the 
Riemann zeta function. 

A functional equation can only be invoked for the semiclassical quantization of 
bound and ergodic systems. In this case the required number of periodic orbits can be 
estimated from the condition s max ~ ng(w) ||, H, which differs by a factor of 4 from 
the required signal length (R3T) for harmonic inversion. Periodic orbit quantization by 
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harmonic inversion will be of particular advantage in situations where special properties 
such as a functional equation cannot be invoked, e.g., for bound systems with non- 
ergodic, i.e., regular or mixed classical dynamics, and for open (scattering) systems. 



4-3. A remark on the Riemann hypotheses 

We conclude this Section with a remark on the famous Riemann hypotheses mentioned 
in Section 2.1. Applying our method of harmonic inversion to the signal C a (s) (Eq. 
18) the Riemann hypotheses for the zeros of the zeta function, £(z = ~ — iw^) = 0, 
is directly related to an equivalent statement for the eigenvalues e~ tTWk of the operator 
e~ ll ~ n , i.e., the generalized eigenvalue problem ( pE|) . Speculations that the operator fl can 
be regarded as the Hamiltonian of a quantum mechanical system have been presented 
by Berry ||. Unfortunately, Cl is not known as it is only defined implicitly by way of 
its matrix representation (28), which is a linear functional of the signal (18). However, 
the very fact that the Riemann zeros are obtained as eigenvalues of some matrix with 
analytically known coefficients is already intriguing. 



5. Periodic orbit quantization 

As mentioned in Section 2 the basic equation (|9|) used for the calculation of Riemann 
zeros has the same mathematical form as Gutzwiller's semiclassical trace formula. Both 
series, Eq. |S] and the periodic orbit sum flHf), suffer from similar convergence problems 
in that they are absolutely convergent only in the complex half-plane outside the region 
where the Riemann zeros, or quantum eigenvalues, respectively, are located. As a 
consequence, in a direct summation of periodic orbit contributions smoothing techniques 
must be applied resulting in low resolution spectra for the density of states [|7], |30f1 . 
To extract individual eigenstates the semiclassical trace formula has to be analytically 
continued to the region of the quantum poles. Here dynamical zeta functions have 
turned out to be of particular interest. 

For bound and ergodic systems one technique is to apply an approximate functional 
equation and generalize the Riemann-Siegel formula (14) to dynamical zeta functions 
0, ||, ||. The Riemann-Siegel look-alike formula has been applied, e.g., for the 



semiclassical quantization of the hyperbola billiard |p9|| . For bound ergodic systems 
alternative semiclassical quantization conditions based on a semiclassical representation 
of the spectral staircase Af(E) = I] n 6(i? — E n ) J7], § and derived from a quantum 
version of a classical Poincare map have also been discussed. 

These quantization techniques cannot be applied to open systems. However, if a 
symbolic dynamics for the system exists, i.e., if the periodic orbits can be classified with 
the help of a complete symbolic code, the dynamical zeta function, given as an infinite 
Euler product over entries from classical periodic orbits can be expanded in terms of 
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the cycle length of the orbits |10|, [TTJ] . The cycle expansion series is rapidly convergent 
if the contributions of long orbits are approximately shadowed by contributions of short 
orbits. The cycle expansion technique has been applied, e.g., to the three disk scattering 
system jn| |12|, jl3fl , the three body Coulomb system ]3T|, ^2|, and to the hydrogen 
atom in a magnetic field A combination of the cycle-expansion method with 



a functional equation has been applied to bound systems in [[H], However, the 
existence of a complete symbolic dynamics is more the exception than the rule, and the 
cycle expansion cannot be applied, in particular for systems with mixed regular-chaotic 
classical dynamics. 

In this Section we apply the same technique that we used for the calculation 
of Riemann zeros, to the calculation of semiclassical eigenenergies and resonances 
of physical systems by harmonic inversion of Gutzwiller's periodic orbit sum for the 
propagator. The method only requires the knowledge of all orbits up to a sufficiently 
long but finite period and does not rely on either an approximate semiclassical functional 
equation, nor does it depend on the existence of a symbolic code for the orbits. The 
method will therefore allow the investigation of a large variety of systems with an 
underlying chaotic, mixed, or even regular classical dynamics. The derivation of an 
expression for the recurrence function to be harmonically inverted is analogous to that 
in Section 2.2. 

5.1. Semiclassical density of states 

Following Gutzwiller |], |2| the semiclassical response function for chaotic systems is 
given by 

g sc (E)=g s c (E)+Y,A po e^° , (37) 

po 

where gQ C (E) is a smooth function and the S po and A po are the classical actions and 
weights (including phase information given by the Maslov index) of periodic orbit 
contributions. Eq. ( p7[) is also valid for integrable and near-integrable |36], [37] 



systems but with different expressions for the amplitudes A po . It should also be possible 
to include complex "ghost" orbits pS), |39| and uniform semiclassical approximations 
p0| , |4l|j close to bifurcations of periodic orbits in the semiclassical response function 
(|57|). The eigenenergies and resonances are the poles of the response function but, 
unfortunately, its semiclassical approximation ([37]) does not converge in the region of 
the poles, whence the problem is the analytic continuation of g sc (E) to this region. 

In the following we make the (weak) assumption that the classical system has 
a scaling property, i.e., the shape of periodic orbits does not depend on the scaling 
parameter, w, and the classical action scales as 

Spo = ws po . (38) 
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Examples of scaling systems are billiards JR], |4"2]| , Hamiltonians with homogeneous 
potentials [|3, Coulomb systems |B3|, or the hydrogen atom in external magnetic and 



electric fields []45|, |33|. Eq. [38] can even be applied for non-scaling, e.g., molecular systems 
if a generalized scaling parameter w = h~g is introduced as a new dynamical variable 
P5fl . Quantization yields bound states or resonances, u>fc, for the scaling parameter. In 
scaling systems the semiclassical response function g sc (w) can be Fourier transformed 
easily to obtain the semiclassical trace of the propagator 

1 r+oo 

C sc (s) = — / g sc (w)e-^dw = £ A po 5 (s - s po ) . (39) 
Air J-oo po 



The signal C sc (s) has 5-peaks at the positions of the classical periods (scaled actions) 
s = s po of periodic orbits and with peak heights (recurrence strengths) A po , i.e., C sc (s) is 
Gutzwiller's periodic orbit recurrence function. Consider now the quantum mechanical 
counterparts of g sc (w) and C sc (w) taken as the sums over the poles Wk of the Green's 
function, 

g qm W = E — djL —~ . (40) 

Y w - w k + xe 

1 f + OO 

(7 qm (s) = — / f(«))e" ,s, "(liD = -iy4e""" feS , (41) 

27T J-oo 7 

with effc being the multiplicities of resonances, i.e., = 1 for non-degenerate states. In 
analogy with the calculation of Riemann zeros from Eq. (18) the frequencies, Wk, and 
amplitudes, dj., can now be extracted by harmonic inversion of the signal C sc (s) after 
convoluting it with a Gaussian function, i.e., 

C s a c (s) = -;L- £ A P J S - S ^ 2 ^ 2 . (42) 
V27r<7 po 

By adjusting C^ c (s) to the functional form of Eq. [41], the frequencies, Wk, can be 
interpreted as the semiclassical approximation to the poles of the Green's function in 
(flOf). Note that the harmonic inversion method described in Section 3 allows studying 
signals with complex frequencies Wk as well. For open systems the complex frequencies 
can be interpreted as semiclassical resonances. Note also that the Wk in general differ 
from the exact quantum eigenvalues because Gutzwiller's trace formula ( ]3~^) is an 
approximation, correct only to first order in %. Therefore the diagonalization of small 
matrices in (Bfj) does not imply that the results of periodic orbit quantization are more 
"quantum" in any sense than those obtained, e.g., from a cycle expansion [0, [llf . The 
eigenvalues are solutions of non-linear equations and the diagonalization is equivalent 
to the search of zeros of the dynamical zeta function in the cycle expansion technique. 
Numerical calculation of the zeros is also a non-linear problem and, in contrast to the 
matrix diagonalization, might encounter a problem of missing roots. 
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5.2. Semiclassical matrix elements 

The procedure described above can be generalized in a straightforward manner to the 
calculation of semiclassical diagonal matrix elements {ip^M^k) of a smooth Hermitian 



operator A. In this case we start from the quantum mechanical trace formula [47 



,7M=trG+i = E ^ fc|A| y. , (43) 

which has the same functional form as ([40]), but with dk = (ipklAfyk) instead of <4 = 1. 
For the quantum response function g ( ^ a (w) (Eq. [43]) a semiclassical approximation has 
been derived in j47j, which has the same form as Gutzwiller's trace formula fl3~7|) but 
with amplitudes 

A po = -i Ap€ (44) 
v /|det(M po -/)| 

where M po is the monodromy matrix and /i po the Maslov index of the periodic orbit, 
and 

A p = f Sp A(q(s),p(s))ds (45) 



o 



is the classical average of the observable A over one period S p of the primitive periodic 
orbit. Note that q(s) and p(s) are functions of the classical action instead of time for 
scaling systems [fS|. Gutzwiller's trace formula for the density of states is obtained 



with A being the identity operator, i.e., A p = S p . When the semiclassical signal C sc (s) 
(Eq. ^9[) with amplitudes A po given by Eqs. ^ and f|5] is analyzed with the method 
of harmonic inversion the frequencies and amplitudes obtained are the semiclassical 
approximations to the eigenvalues Wk and matrix elements dk = (ipk\A\ipk), respectively. 

6. The three disk scattering system 

Let us consider a billiard system consisting of three identical hard disks with unit 
radii, R — 1, displaced from each other by the same distance d. This simple, albeit 
nontrivial, scattering system has served as a model for periodic orbit quantization in 
many investigations in recent years E{I [10], [L2l 0|. After symmetry reduction the 



periodic orbits can be classified by a binary symbolic code fllQfl . For d > 2.1 there is 
a one-to-one identity between the periodic orbits and the symbolic code, whereas for 
d < 2.1 pruning of orbits sets in. For d = 6 semiclassical resonances were calculated 
by application of the cycle expansion technique including all periodic orbits up to cycle 
length n — 13 |Tj|. In order to demonstrate the usefulness of the harmonic inversion 



technique we first apply it to the case R : d = 1 : 6 studied before. Note that the 
ratio corresponds to the very favorable regime for the cycle expansion (see below). In 
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billiards the scaled action s is given by the length L of orbits (s = L) and the quantized 
parameter is the absolute value of the wave vector k = |k| = y / 2mE/h. Fig. 3a shows the 
periodic orbit recurrence function, i.e., the trace of the semiclassical propagator C SC (L). 
The groups with oscillating sign belong to periodic orbits with adjacent cycle lengths. 
To obtain a smooth function on an equidistant grid, which is required for the harmonic 
inversion method, the ^-functions in ( p9|) have been convoluted with a Gaussian function 
of width a = 0.0015. As explained in Section 2 this does not change the underlying 
spectrum. The results of the harmonic inversion analysis of this signal are presented in 
Fig. 3b and Table 3. The crosses in Fig. 3b represent semiclassical poles, for which the 
amplitudes dj. are very close to 1, mostly within one percent. Because the amplitudes 
converge much slower than the frequencies these resonance positions can be assumed to 
be very accurate within the semiclassical approximation. In fact, a perfect agreement 
to many significant figures is achieved for these poles with the results obtained by cycle 



expansion [|T3|| , and this agreement confirms that the results in Fig. 3b and Table 3 
are the true semiclassical resonances, i.e., deviations from the exact quantum poles are 
solely due to the semiclassical approximation in Gutzwiller's trace formula. For some 
broad resonances marked by diamonds in Fig. 3b and Table 3 the c4 deviate strongly 
from 1, within 5 to maximal 50 percent. It is not clear whether these strong deviations 
are due to numerical effects, such as convergence problems caused by too short a signal, 
or if they are a direct consequence of the semiclassical approximation. Of course, in the 
exact expression ( flip all multiplicities dk are 1, but there is no proof that this is still true 
within the semiclassical approximation. However, for the lowest k eigenvalues (see, e.g., 
the first three resonances in Table 3), where the agreement with the exact resonance 
energies is worst |TJ[ d^ = 1 still holds, indicating that there is no ^-dependence for the 
multiplicities. 

The cycle expansion technique is based on the idea that the contributions of long 
periodic orbits are shadowed by those of short orbits. For the three disk scattering 
system this is ideally fulfilled in the limit of a large ratio d/R S> 2. This is no longer true 
for short distances between the disks and hence the convergence of the conventional cycle 
expansion becomes rather slow As a second example of periodic orbit quantization 
by harmonic inversion we therefore study the three disk scattering system with a short 
distance ratio d/R = 2.5, and thus in a situation where the assumption of the cycle 
expansion that contributions of long orbits are shadowed by short orbits is no longer a 
good approximation. The results are presented in Fig. 4 and Table 4. For large L groups 
of orbits with the same cycle length of the symbolic code strongly overlap and cannot 
be recognized in Fig. 4a. The signal is obtained from periodic orbits with cycle length 
n < 13. Note that only 356 periodic orbits with L < 7.5 are included in the signal (Fig. 
4a) whereas the complete set of orbits with cycle length n < 13 consists of 1377 orbits. 
The resonances obtained by harmonic inversion of the semiclassical recurrence function 
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are in good agreement with results of the cycle expansion fl50 |. 

We finally remark that when a symbolic dynamics exists and contributions of long 
orbits are shadowed by short orbits the cycle expansion techniques are very efficient and, 
e.g., for the three disk scattering system with d > 6R a large number of resonances can 
be obtained from just 8 periodic orbits with cycle length n < 4. Again, we do not aim at 
competing with the efficiency of the cycle expansion method in such ideal situations. As 
already mentioned in Section 4 the advantage of periodic orbit quantization by harmonic 
inversion is its general applicability: it does not depend on the existence of a symbolic 
code and the shadowing of orbits. We have extracted the semiclassical resonances of the 
three disk scattering problem directly from the periodic orbits with exactly the same 
method as the Riemann zeros in Section 4 from the primes as 'periodic orbits', for which 
no symbolic dynamics exists. 



7. Conclusion 



We have introduced harmonic inversion as a new and general tool for semiclassical 
periodic orbit quantization. The method requires the complete set of periodic orbits 
up to a given maximum period as input but does not depend on special properties 
of the orbits, as, e.g., the existence of a symbolic code or a functional equation. We 
have demonstrated the wide applicability of the method by applying it to two systems 
with completely different properties, namely the zeros of the Riemann zeta function and 
the three disk scattering problem. Both systems have been treated before by efficient 
methods, which, however, are restricted to bound ergodic systems or systems with a 
complete symbolic dynamics. The harmonic inversion technique allows to solve both 
problems with one and the same method. Therefore the method can also serve as a tool 
for, e.g., the semiclassical quantization of systems with mixed regular-chaotic classical 
dynamics, which still is a challenging and unsolved problem. The signal C sc (s) can 
be composed as the sum of a signal related to the irregular part of the classical phase 
space with periodic orbit amplitudes given by Gutzwiller's trace formula and a signal 
related to stable |35| or nearly integrable 0] torus structures. It should also be possible 
to include, e.g., creeping orbits |HJ, ghost orbit contributions 0, |39], pijfl , and higher 
order h corrections |HJ into the signal C sc (s), which can then be inverted to reveal 
the semiclassical poles. The method can even be used for a semiclassical periodic orbit 
quantization of systems with non-homogeneous potentials such as the potential surfaces 
of molecules when a generalized scaling technique |K| is applied. 
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Figure 1. "Recurrence" function — iC a {s) for the Ricmann zeros which has been 
analyzed by harmonic inversion, (a) Range < s < 6, (b) short range around s = 13. 
The ((-functions have been convoluted by a Gaussian with width a = 0.0003. Dashed 
line: Smooth background C(s) = ie s l 2 resulting from the pole of the zeta function. 
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Figure 2. Estimated number of converged zeros of the Riemann zeta function, which 
can be obtained by harmonic inversion for given number of primes p. 
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Figure 3. Three disk scattering system ( A\ subspace) with R = 1, d = 6. (a) Periodic 
orbit recurrence function, C(L). The signal has been convoluted with a Gaussian of 
width a — 0.0015. (b) Semiclassical resonances. The resonance positions marked by 
diamonds might be less accurate (see text). 
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Figure 4. Three disk scattering system (A\ subspace) with R = 1, d = 2.5. (a) 
Periodic orbit recurrence function, C(L). The signal has been convoluted with a 
Gaussian of width a = 0.0003. (b) Semiclassical resonances. The resonance positions 
marked by diamonds might be less accurate (see text). 



Table 1. Non-trivial zeros Wk, multiplicities dk, and error estimate e for the Riemann 
zeta function. 
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Table 2. Trivial zeros and pole of the Riemann zeta function. 
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Table 3. Semiclassical resonances, multiplicities, and error estimates for the three 
disk scattering problem {A\ sub-space) with R = 1, d = 6. The marked resonances are 
plotted as diamonds in Fig. 3b. 
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Table 4. Scmiclassical resonances, multiplicities, and error estimates for the three 
disk scattering problem (A\ subspace) with R = 1, d = 2.5 obtained from the signal 
C(L) with L < 7.5 in Fig. 4a. The marked resonances are plotted as diamonds in Fig. 
4b. 
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